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Abstract 



The launching process of a magnetically driven outflow from an accretion disk is investigated in a local, shearing box model which 
allows a study of the feedback between accretion and angular momentum loss. The mass-flux instability found in previous linear 
analyses of this problem is recovered in a series of 2D (axisymmetric) simulations in the MRI-stable (high magnetic field strength) 
regime. At low field strengths that are still sufficient to suppress MRI, the instability develops on a short radial length scale and 
saturates at a modest amplitude. At high field strengths, a long- wavelength "clump" instability of large amplitude is observed, with 
growth times of a few orbits. As speculated before, the unstable connection between disk and outflow may be relevant for the time 
dependence observed in jet-producing disks. The success of the simulations is due in a large part to the implementation of an effective 
wave-transmitting upper boundary condition. 
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1. Introduction 

The occurrence of strong, collimated outflows in association 
with the accretion of compact astronomical objects is a common 
phenomenon. Examples are jets from young stellar objects, in- 
teracting binaries, active galactic nuclei and possibly the central 
engines of gamma ray bursts. Although they are often regarded 
as separate entities, accretion disks and jets are physically de- 
pendent on one another. The mass loading of the jet depends on 
properties of the disk and its immediate atmosphere, and the out- 
flow feeds back on the disk by carrying off material and angular 
momentum. The physics of the transition from accretion flows 
to outflows is therefore a key element of our understanding of 
both disks and jets. 

An inflow of gas in a rotating disk necessitates the pres- 
ence of a mechanism for angular momentum transport. It is well 
known that weak magnetic fields are responsible for magne- 
torotational instabilities (MRI; Balb us & Hawley||1991| [Balbus 
2003). The turbulence driven by MRI leads to internal stresses 
that redistribute angular momentum in a similar way as would do 
a high viscosity, thus driving accretion. MRI in disks have been 
the subject of numerous numerical studies, including local sim- 
ulations in shearing boxes as well as global simulations (e.g., 
Hawley & Krolik|[2002l |Turner et aLp002| |Fromang & Stone| 



2009) . 



Accretion disks may function without MRI, or other forms of 
viscosity, if they lose angular momentum through outflows. Such 
outflows are generated r eadily in disks threaded by large-scale 
poloidal magnetic fields ([Bisnovaty i-Kogan & Ruzmaikin|1976j 
Blandfor d|1976| ). The material is accelerated outwards along the 
poloidal field by centrifugal forces. The inclination of the field 
with respect to the disk determines the efficiency of the pro- 
cess. It must be sufficiently small, <60° in the cold wind model 
of Blandford & Payne| ( |1982 ), for the "magnetic slingshot" to 
work. In an optically thick disk with an isothermal atmosphere 
and a mean field that is sufficiently strong to suppress MRI, the 
strongest winds are expected for inclinations of ^45° and the 



mass loss rate decreases with increasing field strength; in turbu- 
lent disks with weaker fields, the mass loss rate is expected to 
increase with the field strength (Ogil vie & Livio|2001| ). 

Due to the way in which the inclination of the poloidal mag- 
netic field is coupled with the strength of the outflow, the connec- 
tion between accretion and magnetically driven outflows might 
be inherently unstable. A possible instability scenario is the fol- 
lowing: an increase of the speed of radial inflow increases the 
inclination of the field towards the midplane, which increases 
the mass flux and the loss of angular momentum and thus leads 
to an even faster inflow. The existence of instabilities in disks 
that lose angular momentum through a magnetically driven wind 
was first predicted by |Lubow et al.| ( |1994 ), disputed by Konigl 
& Wardl e](|1996|), an d co nfirmed by the perturbation analyses of 
Cao & Spruit| ( |2002| ) and |Campbell| ( |2059l ). 

The perspective on the physics of the disk-wind connection 
taken here starts with the stationary, one-dimensional flow prob- 
lem of |Qgilvie & Livio] ( |2001| ), which yielded the dependence 
of the outflow rate on the strength and inclination of the mag- 
netic field. A logical step towards a time-dependent view are lin- 
ear stability analyses of this stationary problem, such as the ax- 
isymmetric, short- wavelength problem studied in |Cao & Spruit 
( 2002 ). The logical step taken here is a two-dimensional, finite 
amplitude study of the same problem by numerical simulations. 

Borrowing from early studies of MRI, the problem is kept lo- 
cal in the radial direction by using a Cartesian shearing box. The 
simulations can thus be seen as an extension of 2D MRI simula- 
tions with a net magnetic flux threading the box. The extension 
consists of including an outflow through the upper boundary, 
and a finite asymptotic inclination of the magnetic field. In this 
way, it provides a natural connection between MRI and wind- 
launching studies like those of Ogilvie & Livio. 

The model is described in Sect. [2] Details of the numerical 
realization are given in Sect. [3] The results of the simulations are 
presented in Sect. |4| followed by a summary and discussion in 
Sect.E 
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2. The model 

We consider an axisymmetric accretion disk threaded by a 
poloidal magnetic field that is bent away from the rotation axis. 
The accreting matter, which is initially orbiting the central ob- 
ject in (near) radial force equilibrium, loses angular momentum 
through a centrifugally accelerated wind, thus enabling an ac- 
cretion inflow. We ignore the possibility of angular momentum 
transport by viscosity or MRI, and, with one exception, assume 
no explicit magnetic diffusivity. 

The atmosphere of the disk is assumed to have a similar tem- 
perature as the disk itself. Both density and plasma-/? decrease 
strongly with height. Inside the disk, J3 is of order unity or larger, 
and the frozen-in magnetic field is dragged along with the rotat- 
ing plasma. Above the disk, the magnetic field is strong enough 
to enforce an approximately constant angular velocity along the 
field lines. The material trapped on the magnetic field is flung 
outwards by centrifugal forces. Around the point where the flow 
reaches Alfven speed, the inertia of the gas dominates again and 
causes a strong winding of the magnetic field into a predomi- 
nantly azimuthal field. 

The simulations are done in a Cartesian, periodic shearing 
box (see, e.g., the MRI simulations by Hawley et al.|1995| ). The 
focus is thus on local processes, happening on a radial scale 
of a few times the thickness of the disk and smaller. The goal 
is to study the nonlinear development of short-wave processes 
such as the instabilities predicted by linear analyses of the wind- 
launching problem. A major advantage of the shearing box is 
that the calculations can be done without adding an explicit mag- 
netic diffusivity. With (quasi-)periodic boundary conditions in 
the radial direction, magnetic flux does not pile up by advection, 
as would happen in a global model. In effect, the shearing box 
model thus elegantly incorporates a separation of time scales: it 
includes processes that scale with the orbital period while leav- 
ing out long-term processes such as the pile-up of magnetic flux. 

Since a wind is to be launched, a realistic upper boundary 
that allows mass outflow is needed. At low mass fluxes, the mag- 
netic field configuration above the disk is affected only little by 
the presence of the outflow. It is therefore natural to approximate 
the poloidal field at the upper boundary as a potential field. The 
implementation of this boundary turns out to be critical for the 
success of the model. 

The mass loss from the box causes conditions to change 
slowly with time. In order to separate such secular changes from 
the processes under study, we add a mass sourc e clos e to the 
midplane, as in previous local models ( |Ogirvie||2012| ). As an 
approximation, this should be valid under conditions where the 
time scale for mass loss is long compared with the accretion time 
scale. 



2.1. Magnetic support against gravity 

The poloidal field in this model is not force-free. Near the mid- 
plane, the bending of the field generates a curvature force which 
counteracts gravity. The midplane equilibrium is thus deter- 
mined by inward gravity, outward centrifugal force and outward 
magnetic curvature force. Let 
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be an estimate for the relative importance of the curvature force 
at the midplane if the rotation is Keplerian. Taking the curvature 
radius r c to be of the order of the disk's scale height H and es- 
timating Q 2 with c 2 /H 2 = p/H 2 p (thin disk approximation), we 



find 6 ~ S/J3. For J3 > 1 and a small aspect ratio 6 = H/R, the 
curvature force has only minor significance. It perturbs the radial 
equilibrium of a Keplerian disk only slightly. This difference is 
important, however, for the launc hing of the magnetically pow- 
ered wind ( Ogilvi e & Livio|1998] ). 



3. Methods 



We use the shearing box approach (e.g., |Hawley et al.|1995| ) with 
axisymmetry to solve the MHD equations in a Cartesian box that 
rotates with Keplerian angular velocity Qo = Q$ z at some dis- 
tance Ro from the axis of rotation. Compared to the dimensions 
of the box, Ro is assumed to be large enough that the equations 
to be solved become independent of it. 

The majority of simulations was done assuming ideal MHD, 
for which the induction equation has the usual form d t B = V x 
(v x B). The momentum equation is 

dv _ Vp (VxB)xB 

dt p 4np 

- 2A) x v + Q l oOxe x - ze z ) - -xQ M, (2) 

where x = R - Rq is the radial coordinate and z is the height (ver- 
tical) coordinate, z = corresponding to the midplane. The first 
two terms after the Lorentz force are the result of a Taylor ex- 
pansion of the centrifugal force RQ^e R and the gravitational ac- 
celeration -Q 2 R^r~ 2 e r , assuming that \x\ , |z| <^ ^o- The last term 
accounts for the momentum of material which is added through 
a mass source term (described below). We also calculated diffu- 
sive cases in which the induction and momentum equations are 
extended by r]V 2 B and vV 2 v, respectively. 

To prevent the box from being drained by the outflow, we 
introduce a mass source term on the right-hand side of the con- 
tinuity equation: 



— + V • (pv) = f 

dt r M 



M 



(3) 



where / = 1 for \z\ < 0.5, / = cos[tt(|z| - 0.5)] for 0.5 < \z\ < 1 
and / = for \z\ > 1. Integrated in z, the material introduced 
in the time span tm amounts to Zm/^o = 61% of the surface 
density in the initial state. It is given the initial temperature and 
Keplerian momentum (terms with M in Eqs.|2]and|4]). It damps 
horizontal and vertical motions and stabilizes the temperature. 

We adopt an ideal gas equation of state with y = 5/3 and 
evolve the equation for the internal energy density e = p/(y-l), 



de e(t = 0) 

+ V • (ev) = -pV • v + K + M—j- — — ^, 
dt pit ~ 0) 



(4) 



along with the corresponding equation for conservation of total 
energy. In places where the latter yields unphysical results ("neg- 
ative pressures" due to discretization errors which occur occa- 
sionally in very low-/? regions), we use the value obtained by 
evolving the internal energy instead of the total energy. The last 
term in Eq. ([4]) accounts for the material which is added through 
the mass source term. Since a calculation can become unfeasible 
if regions with extremely low density develop, we add another 
artificial source term 

T{t = 0) - 7X0 e(t = 0) 



K = 



Tit = 0) 



(5) 



to the energy equation. This intervention helps to avoid low- 
density cavities by relaxing the temperature to that of the initial 
state on an appropriate time scale tr. In nature, radiative losses 
would prevent excessive temperatures. 
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3.1. Boundary conditions 

We assume reflective symmetry at the midplane. The computa- 
tional domain is -L x < x < L x in the horizontal direction and 
< z < L z in the vertical direction. The bottom boundary corre- 
sponds to the accretion disk's midplane. There, we use reflective 
boundary conditions: d(p,T 9 V[ Xt y],\v z \,\B[ Xt y]\)/dz = 0, the signs 
of v z and B[ x ^ are reversed and B z follows from solenoidality. 

The horizontal boundaries are strictly periodic for all quan- 
tities except the azimuthal velocity, for which 



v y (x, z) = v y (x ± 2L X , z) ± -L x Qq 
is applied in the ghost cells of the left (right) boundary. 



(6) 



3.1.1. Top boundaries 

The top boundary conditions (z > L z ) must account for the ef- 
fects of a global magnetic field outside the scope of the local 
simulation and allow for an unhindered outflow of material. We 
tried different ways of implementing these constraints and found 
many to be unfit: zero-slope conditions in the vertical direc- 
tion create numerical instabilities and rigorously imposing the 
poloidal field inclination angle causes significant reflections, in 
some cases strongly affecting the dynamics in the simulated do- 
main. 

The conditions described below are constructed in view of 
the limiting case far above the midplane. There, J3 «: 1, which 
suggests that the magnetic field is force-free, and the characteris- 
tics of MHD waves along the field are outwards, which suggests 
the use of extrapolation along the magnetic field. A convenient 
assumption then is that the poloidal magnetic field is a potential 
on^] An inclination angle is imposed by choosing the value of 
the mean field. This gives the field more freedom than a strict 
imposition of the inclination angle. These boundary conditions 
work very well even in cases where the boundary is not far in the 
<$c 1 domain, still inside the Alfven surface. They are numeri- 
cally stable and cause only minimal reflective artifacts. 

The conditions are implemented as follows. At each time 
step, we determine the potential field that matches with B z at 
z = L z and satisfies [V x B] y = 0. The solution is found 
by means of a discrete Fourier transform in x-direction, using 
B x (k) = -iB z (k) for the complex Fourier coefficients of B x (B z 
being the coefficients of B z ) for k ^ 0. The mean field (corre- 
sponding to k - 0) is chosen such that at infinity the poloidal 
field is inclined by an angle i with respect to the horizontal. We 
use zero- slope conditions along the magnetic field: with s mea- 
suring distance along a poloidal field line, dip, V[ Xf y fZ ],B y )/ds = 
to first order (i.e., the values are constant along a field line). In 
addition, we dismiss inertial and gravitational forces parallel to 
the magnetic field in the uppermost layer of cells next to the up- 
per boundary. The temperature is kept fixed at the top boundary, 
dT/dt = 0. 

3.2. Initiai conditions 

The gas is initially isothermal, with both pressure and density 
being oc exp[-z 2 / (2/^)], where /q is the unit length (i.e., /q is the 
scale height in the initial state; we shall denote the actual, time 
dependent scale height of the disk with H). The initial magnetic 



1 It should be borne in mind, though, that this may not be a good 
approximation in cases where the azimuthal field exerts strong Lorentz 
forces perpendicular to the poloidal field. 
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Figure 1. Visualization of the wind properties listed in Table 
Each "flower" represents a simulation. The radii of the "petals" 
are proportional to the represented value. All values are normal- 
ized by the respective value in run /3U45. Strong clumps develop 
in the cases to the left and above the dotted line (cf. Figs. [4] [5] 
and [6]). Note that the j3 in the simulation identifiers refers to the 
initial state and measures the absolute field strength, whereas the 
J3 represented by the gray petals refers to the actual value in the 
evolved state (see Table [T] and text). 



field is homogeneous and inclined by an angle i with respect to 
the midplane. Its strength is determined by the simulation param- 
eter J3{, which represents the ratio of the gas-to-magnetic pres- 
sure at the midplane in the initial state. The initial velocity is 
KepleriarJ^] at the midplane, = -3/2x£2 , and constant along 
the magnetic field lines. 



3.3. MHD code 

The simulations were performed with the Eulerian MHD code 
of |Obe rgaulinger ( 2008 ). It is based on a flux-conservative 
finite-volume formulation of the MHD equations and a con- 
straint transport scheme that maintains a divergence-free mag- 
netic field ( [Evans & Hawley 19 88]). Using high-resolution shock 
capturing methods (e.g., LeVeque| [l992), it allows a choice of 
various optional high-order reconstruction algorithms and ap- 
proximate Riemann s olvers based on the multi-stage method 
( |Toro & Titarev|2006] ). The simulations presented here were per- 
formed using a fifth order monotonicity -preserving reconstruc- 
tion scheme (Suresh & Hu ynh|1997| ), the HLL Riemann solver 



( |Harten|1983 



, and third order Runge-Kutta time stepping. 



4. Results 



We performed several simulations with varying physical and nu- 
merical parameters. Long-lived runs, in which the disk evolution 
could be followed over many orbits past an initial transient, are 
listed in Table [T] The time scale tk for the temperature con- 
trol was always taken to be 0.1 orbits, so that large deviations 
from the initial isothermal state are rare; e.g., in run /3k45, only 
about 4% of the volume deviates from the initial temperature by 
more than 10%. Unless noted otherwise, the time scale tm for 



2 Taking magnetic support against gravity (Sect. into account 
would change v y by a constant amount -eR Qo/2. We opted not to in- 
clude such a correction and let the system find an equilibrium by itself 
instead. 
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Table 1. List of simulations and mean properties in the evolved state 



Run a 


Domain 


Resolution 


a - «p)) 


((E)) 


((H)) 




((pV.AVy)) 


((-B z B y /4n)) 


((m-m)) 


Pol 


" «Pmag» 


K>] 






[«s»«f0>^] 




[m m Q ] 


pi i45 


4x2 


256 x 128 


2.76+8% 


0.718+3% 


0.651+3% 


0.0677+29% 


0.0897 ±65 % 


1.41+13% 


3.07 ± i7% 


£1*45 h 


4x2 


512x256 


2.43 ±8 % 


0.642+2% 


0.668+6% 


0.0756 ±28 % 


0.108+82% 


1-15+25% 


2.56+60% 


pi l45 hv 


4x2 


512x256 


2.55+4% 


0.695+4% 


0.668+2% 


0.0687+8% 


0.0807+10% 


1.27 ±7 % 


2.71+6% 


J31i45h7] 


4x2 


512x256 


3.19~ 6 % 


0.789 ± i% 


0.639~ 4 % 


0.0615+20% 


0.0865 ±76 % 


1.61+8% 


3.48+i9% 


pi *45 b 


6x3 


384 x 192 


2.56+20% 


0.624+i % 


0.661+15% 


0.0744+29% 


0.111 ±72 % 


1.31+27% 


2.97+34% 


p\i31 


4x2 


256 x 128 


3.76+i8% 


0.586 ±3 % 


0.557~ 13 % 


0.0831+39% 


0.177 ±55 % 


1.93+20% 


4.27+40% 


pi i40 


4x2 


256 x 128 


3.27+i6% 


0.626+2% 


0.583+9% 


0.0775+23% 


0.143+42% 


1.70 ±2 i% 


3.79 ±27 % 


pi i50 


4x2 


256 x 128 


1.54+i4% 


0.762+io% 


0.723+9% 


0.0677 ±83 % 


0.118+174% 


0.617+30% 


1-51+59% 


pii60 


4x2 


256 x 128 


1.89+i4% 


1.15+n% 


0.944+5% 


0.0428 ±60 % 


0.0456 ± no% 


0.190+59% 


0.494 ± i83% 


P0.5 i45 


4x2 


256 x 128 


1.69*35% 


0.909 ± i 6 % 


0.627+20% 


0.0529 ± i35% 


0.095 l ±3 i7% 


1.21 + 104% 


2.64+i4o% 


P2i45 


4x2 


256 x 128 


3.46 ±6 % 


0.486+2% 


0.698+3% 


0.100+15% 


0.147 ±3 i% 


1.05 ± n% 


2.46 ± i 5 % 


P% i45 


4x2 


256 x 128 


7.12+i% 


0.286+o% 


0.793 ± i% 


0.170+2% 


0.219 ±4 % 


0.447+2% 


1.34 ± i% 


/?10*45 


4x2 


256 x 128 


7.66+9% 


0.270+2% 


0.806+3% 


0.180+8% 


0.227 ± i 3 % 


0.480+7% 


1.42 ±8 % 



Notes. The quantities in columns 4-10 are average values in the evolved state (for a visualization, see Fig. [I}. Double angle brackets ((. . .» denote 
a combined temporal and horizontal mean. Temporal fluctuations are rounded to full percents. (a) The label contains the initial value of pi p mag at 
the midplane (denoted Pi in the text) and the asymptotic inclination angle i of the potential field at the outer vertical boundary (0° is horizontal and 
90° is vertical). In addition, "h" stands for high resolution, "v" stands for the use of an explicit shear viscosity, "77" stands for the use of explicit 
magnetic diffusion and "b" stands for a bigger computational box. 
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Figure 2. Horizontal means of the surface density (green line), 
scale height (blue line), mass inflow rate (red line), vertical mass 
flux (cyan line), Reynolds stress (brown line) and Maxwell stress 
(purple line) as functions of time in run filt45. 



the mass source used was 2 orbits (= 2n/Qo = to, used as unit 
time). In a simulation similar to fili45 but without a mass source 
(tm —> 00), the mass inside the box is diminished by 58% dur- 
ing 2 orbits (after which the increased Alfven velocities make it 
impractical to continue the simulation). The material introduced 
with the tm = 2 mass source is roughly of this value (cf. Em/Eq 
in Sect. [3]). The simulations were continued to ~20 orbits, at 
which time the system has passed through the initial transient 
and has entered a (more or less, depending on the parameters) 
quasi-stationary state. 

The two basic simulation parameters are Pi, which deter- 
mines the strength of the magnetic field, and t, which determines 
the asymptotic inclination of the poloidal magnetic field in the 
disk's magnetosphere. Pi is usually smaller than the mean value 
of p at the midplane in the evolved state, which is listed in the 
4th column of Table [T] An exception are disks with very weak 
magnetic fields, see Sect. |4.5| 

Fig. [2] shows the evolution of horizontally averaged quanti- 
ties in run pii45. In general, the simulations start with strong 
epicyclic oscillations in the inflow rate, presumably governed by 



1 .4 
1.2 
1 .0 
0.8 
0.6 
0.4 





I 


P iPo] - 






■ P / 2.7 






■ c s [I0l /t ] 






v c [10 I /t ] 






■ v Ap [30 l /t ] 






■ -B y /B z 












X \ \ 

^^Y^-- ^^^^ 









Figure 3. Vertical stratification of the disk in run /3U45: density 
(green line), plasma-/? (blue line), sound speed (oc square root 
of temperature, red line), slow magnetosonic cusp speed (cyan 
line), poloidal Alfven velocity (brown line) and azimuthal mag- 
netic field (magenta line). The solid and dashed black lines rep- 
resent the vertical velocity and the velocity along the poloidal 
magnetic field, respectively. The curves represent the mean over 
the horizontal direction and time in the evolved state. For com- 
parison, dotted lines represent the isothermal, Gaussian stratifi- 
cation of the initial condition (in arbitrary normalization). 



the interplay of magnetic curvature forces (d ue to the bending 

, the Coriolis 



of the magnetic field at the midplane, see Sect. 2.1 
force and the loss of angular momentum by the magnetically 
powered wind. The disk gradually enters a new equilibrium state 
with a vertical outflow. For further analyses, we consider only 
the times after the initial transient. Depending on the simulation, 
the evolved state is more or less relaxed. In some cases, the os- 
cillations are still strong after many orbits (see fluctuations given 
in percent in Table [T}, with no decreasing trend at the end of the 
simulated time span. 

The green curve in Fig. [2] shows the surface density E = 
J pdz, normalized by its initial value. We consider it a mea- 
sure for the efficiency of the wind, with strong winds yielding 
a smaller E. In the evolved state, the mass loss through the wind 
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t=1 7.93 




Figure 4. Snapshot of the density (black/white) and the field- 
parallel mass flux (green/pink colors) in run /3U45. Selected 
magnetic field lines are drawn in blue. The additional solid and 
dashed lines indicate the scale height and the sonic surface, re- 
spectively. 



is balanced by the mass source, so that E « constant. The blue 
curve shows the scale height H , which we define here as the 
height that contains 68% of the mass (we take into account the 
finite vertical size of the computational box, so that H = Iq in the 
initial state). The cyan curve shows the vertical mass flux near 
the upper boundary, measured at height z = 1.7. In the evolved 
state, its mean value equals the amount of material that is intro- 
duced by the mass source: pv z « Em/tm = 0.0616 EQo. The 
brown and purple lines show the Reynolds stress pv z Av y , where 
Av y denotes the deviation from the Keplerian rotation velocity, 
and the Maxwell stress -B z B y /4n, respectively. The red curve 
shows the mass inflow rate m m - J pv x dz, normalized such that 
it gives, in units of the scale height, the horizontal width of the 
disk material that crosses a surface of constant radius per radian 
rotation. 

Fig. [3] shows (horizontally averaged) vertical profiles of 
various quantities in run J31l45. Density and plasma-/?, repre- 
sented by the green and the blue line, decrease by factors of 
15 and 50, respectively, from the midplane to the upper bound- 
ary. Both quantities drop considerably faster with height than 
in a Gaussian stratification. The mean of the sound speed, rep- 
resented by the red line, is close to the initial value (which is 
yy27r/o/^o) at every height. Above the disk, it approximates the 
slow magnetosonic cusp speed v c , given by v\ - c\v 2 h j{c\ + v^) 
and represented by the cyan line. The Alfven velocity of the 
poloidal field, va p , is represented by the brown line; it is by a 
factor of four larger than the cusp speed at the upper boundary. 
The magenta line represents the azimuthal magnetic field com- 
ponent, normalized by the vertical field (the mean of which is a 
constant). 

Fig. [4] shows snapshots of the density p (in units of the initial 
value at the midplane) and of the mass flux pv\\ in the direc- 
tion of the poloidal magnetic field (v\\ is computed from v z and 
the local field inclination) in run J31i45. The outflow is inhomo- 
geneous, the mass flux varies across different field lines. Dark- 
green stripes indicate the regions where the wind is strongest. 
Along some of the field lines, there is an inflow at low heights 
(pink colors). As angular momentum is lost through the wind, 



an accretion flow drags the magnetic field inward (i.e., in -x- 
direction). A time animation shows that the stripes move inward 
with the magnetic field. On a given field line, the strength of the 
wind varies periodically with a frequency of -1.2 per orbit. 

4.1. Dependence on field strength and inclination 

The temporal means of the quantities discussed above are listed 
in Table[T]for different simulations. A graphical representation of 
these numbers can be found in Fig. [T] Snapshots of simulations 
with different values for the asymptotic poloidal field inclination 
l are collected in Fig. [5] those with different values of J3{ are 
collected in Fig. [6] 

With increasing field inclination i, the strength of the wind, 
as measured by the value of the surface density E or the density- 
normalized mass flux pv z , decreases. The disk is less compact 
(increasing H) and the inflow rate smaller. The disks in runs with 
high inclination, l = 50° and i - 60°, develop a large-scale in- 
stability in the form of a density clump. We discuss the clump in 
the i = 50° case in Sect.l4~4l 

With increasing fi x (decreasing field strength), E decreases, 
hence the wind becomes more effective at removing disk mate- 
rial. The inflow rate m m , however, is smaller. While the Reynolds 
stress increases, the inflow rate and the Maxwell stress, which is 
the dominant agent for the removal of angular momentum, both 
decrease. The wind is more homogeneous across the magnetic 
field in the high-/? cases. 

The mean value of ft at the midplane is usually higher than 
the initial value f3{ as a result of the arbitrary initial transient. For 
equal asymptotic inclinations i, the field inclination as a function 
of height is different in simulations with different field strengths. 
For instance, in run J38i45, the horizontal separation between the 
foot point of a field line at the midplane and the intersection 
point of the same field line at the upper boundary is typically 
larger than in run /3U45, see Figs. [5^ and [6^. That is, the mean 
inclination of the field is higher in run /38l45. In general, the 
inclination is a function of the field line as well as height (or 
distance along the field line). 

4.2. Resolution, diffusion and box size 

A comparison of simulations f3li45 and /3k45h, Figs. [4] and [TJ), 
shows that the inhomogeneities (stripes) in the mass flux are nar- 
rower at higher resolution. The short wavelengths appear to grow 
fastest, with numerical resolution the limiting factor. 

We study the effects of diffusion of magnetic field and mo- 
mentum through a high-resolution simulation (J3li45hr]) with an 
explicit magnetic diffusivity rj = lO -3 ^/^ (= 10" 3 c 2 /Qo with 
CsO = OWpo) 1/2 being the isothermal sound speed in the ini- 
tial state) and another high-resolution simulation (/3k45hv) with 
an explicit shear viscosity v = 10~ 2 QqI 2 . In the case with mag- 
netic diffusivity, Fig.|7]l (J31l45\\t]), the stripes are broader than 
in the low-resolution simulation without diffusivity. In the case 
with viscosity, Fig. [7]: (J3li45hv), the wind inhomogeneities are 
mostly gone but the disk develops a clump. The disk is also more 
affected by instabilities in a simulation with a bigger domain 
(/3U45b), see Fig. [7^. The mean wind properties, however, are 
similar in all these cases, see Fig. [T] We conclude that the am- 
plitude of the instability increases with length scale, its shortest 
length scale is limited by (real or numerical) diffusion, its largest 
length by the box size. 
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R. Moll: Shearing box simulations of accretion disk winds 




(b) J31 *40 

2.0 




Figure 5. Snapshots of the density (black/white) and the field-parallel mass flux (green/pink colors) in the evolved state of simula- 
tions with different values for the asymptotic field inclination: (a) i - 37°, (b) i - 40°, (c) i - 50°, (d) i - 60°. The i - 45° case 
is shown in Fig. [4] The plots include selected magnetic field lines, the scale height (solid horizontal lines) and the sonic surface 
(dashed lines). 



Table 2. Average value of the surface density 2 [2o] in runs with 
different tm 





(fc,*) = (l,50°) 


(fc,0 = (M5°) 


(fc,0 = (M5°) 


0.1 


5.32 


5.18 


3.91 


2 


0.762 


0.718 


0.286 


4 


0.524 


0.481 


0.178 


10 


0.376 


0.259 


0.106 



Considering runs with different field strengths and inclinations, 
2 varies in the same direction for every value of tm- 

With increasing values of tm, the solutions become more un- 
stable, developing the typical stripes and clumps. Figure[8]shows 
this for the standard case. As an increase of tm leads to a de- 
crease of density and gas pressure, this is consistent with the 
observation that the amplitude of the instability increases at low 
P. 



4.3. Replenishment of mass 

In the simulations presented above, the value of tm is always 2. 
This not only facilitates the comparison, but it also makes practi- 
cal sense because it keeps the surface density at values which are 
not very far from the initial state (which is that of a classical thin 
disk). Since the mass source is somewhat artificial, it is reason- 
able not to have it too strong, where "strong" means putting in 
more mass than what would be lost without it in the course of a 
few orbits. On the other hand, simulations with very low density 
are unfeasible for numerical reasons. 

Table [2] shows how the surface density changes for different 
values of the tm- As expected, it increases with decreasing tm- 



4.4. Clumpy disks 

In some simulations, strong horizontal inhomogeneities in the 
density develop: the disk condenses into one or several persistent 
"clumps". A single distinctive clump forms in run J31l50, see 
Fig. [3£. During the formation, the outward mass flux is often 
higher along field lines which are not anchored in the clump, see 
Fig. [9] The clump's rotation is relatively slow (sub-Keplerian). 
In its wake, plasma-/? is high and the field at low heights is highly 
inclined towards the midplane. 

We measure the growth of the clump by considering the in- 
homogeneity of the inclination of the poloidal magnetic field in- 
side the disk, at a height of z = 0.3 « 0A((H)). The result is 
depicted in Fig.flO] The perturbation grows exponentially with a 
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(a) /?2*45 




(b) £0.5*45 

2.0 



(c)/?8*45 

2.0 





(d) £10*45 

2.0 




Figure 6. Snapshots of the density (black/white) and the field-parallel mass flux (green/pink colors) in the evolved state of simu- 
lations with different values for the initial gas-to-magnetic pressure: (a) £i = 0.5, (b) £i = 2, (c) fii = 8, (d) fii = 10. The fii = 1 
case is shown in Fig. [4] The plots include selected magnetic field lines, the scale height (solid horizontal lines) and the sonic surface 
(dashed lines). 



growth time of 2.5 orbits. Saturation is likely influenced by the 
horizontal periodicity of the computational domain. In simula- 
tion /30.5*45, which also develops a single clump, we measure 
an exponential growth time of 4.2 orbits by the same method. 

Fig. [TT] shows the evolution of various quantities, (a) along a 
field line for which the density at low heights decreases during 
clump formation and (b) along a field line for which the density 
increases. The difference in density relates to an overall stronger 
mass flow (green colors) along field line (a). Line (b) develops a 
strong twist B y /B z and a high inclination towards the midplane 
at low heights. Rightwards inclined stripes in the mass flux plot 
indicate outward moving perturbations in the outflow. There are 
also fan- shaped features which indicate perturbations that move 
down from the top and are reflected at the midplane. 

4.5. Weak magnetic fields and MR I 

The upper magnetic boundary conditions, viz., a potential 
poloidal field, are introduced under the assumption of strong 
magnetic fields with J3 <?c 1 . In the initial stratification, J3 - 1 
at z = in run j3li45 and j3 = 1 at z = 2.1 (i.e., slightly above 
the upper boundary) in run f3 10*45, which, of all cases presented 
above, is the simulation with the weakest field. The respective 



values for ((/?)) /((/? m ag)) in the evolved state are z = 0.32 and 
z = 0.44; both are well below the upper boundary. That is, with 
the exception of perhaps an initial transient, the boundary con- 
ditions are physically sound in the cases presented above. 

The orderliness of the magnetic field is destroyed if the sys- 
tem is subject to MRI. As expected, we see such cases in simu- 
lations with very low field strengths. Fig. [T2| shows snapshots of 
two simulations with f3{ = 100 and no mass source (tm — > °o). In 
both runs, the initial azimuthal velocity v y was perturbed by ran- 
dom perturbations with an amplitude of 1 c s o ~ 6.3 /o/fo- m the 
first case (top panel), * = 90°. The initially vertical magnetic field 
is radially stretched by axisymmetric MRI modes (compare, e.g., 
|Stone & Gard iner 2010). The perturbations grow strong within 
roughly an orbit. In the second case (bottom panel), * = 45°. The 
magnetic field becomes very flat at low heights and the density 
develops a strong peak: (p(z = Q))/(p(z = 0.3)) « 12 at t - 0.97. 
The field then reconnects at the midplane. 

Runs with * = 45° and J3 { = 12, 100, 10 20 are all numerically 
unstable (i.e., some values become too extreme for the numeri- 
cal solver) and terminate after -1/2 orbit. Common to all these 
cases is that the magnetic field becomes very flat near the mid- 
plane. Unlike in the low-/? cases, the reflecting condition does 
not make the field vertical near z = 0. For J3{ — > oo (no magnetic 
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(a)/?U45b 

3.0 



(b)£l *45h 





Figure 7. Snapshots of the density (black/white) and the field-parallel mass flux (green/pink colors) in the evolved state of variants 
of run f3li45 shown in Fig. [4] The panels show the effect of (a) a bigger computational box, (b) a doubled resolution, (c) a doubled 
resolution and explicit shear viscosity, (d) a doubled resolution and explicit magnetic diffusivity. The plots include selected magnetic 
field lines, the scale height (solid horizontal lines) and the sonic surface (dashed lines). 



field), as expected, the disk is stationary and does not generate a 
wind. 

To test whether the peculiar behavior at very weak magnetic 
field strengths is tied to the initial state, we continued simulation 
J3li45 with the magnetic field strength reduced by a factor 10. 
A snapshot of this run is shown in Fig. 13 Due to the changed 
magnetic forces, the originally stable system is out of equilib- 
rium and oscillates. The magnetic field is distorted by instabil- 
ities at all heights. After about 1 orbit, an elongated magnetic 
structure develops near the midplane, moving inward fast. The 
mean plasma-/? is smaller than 1 for z > 1.35, i.e., the potential 
condition at the upper boundary is still practical in this simula- 
tion. 



5. Summary and discussion 

We have presented simulations of magnetocentrifugally acceler- 
ated winds in an axisymmetric shearing box. The box contains a 
thin disk with an ordered poloidal magnetic field that is inclined 
away from the rotation axis. Special upper boundary conditions 
were used that allow a study of the coupling between the disk 
and the angular momentum removing outflow. 



For magnetic fields which are strong enough to suppress 
MRI, we find that material is efficiently accelerated away from 
the disk. To be able to follow the disk evolution over many or- 
bits, we replace lost material by a time-independent mass source. 
After several orbits, the wind enters a quasi- steady equilibrium 
with the mass source. The strength of the wind, as measured 
by how efficiently it limits the amount of material in the box 
in spite of a constant resupply, increases with smaller magnetic 
field strengths or smaller asymptotic inclinations of the poloidal 
magnetic field. Angular momentum is lost mainly through mag- 
netic torques (Maxwell stresses). The disks develop a radial in- 
flow with sub-Keplerian rotation velocities near the midplane. 

In general, the outflows are not homogeneous and the mag- 
netic field lines are not uniform. All cases develop instabilities in 
the mass flux on the shortest length scales that are numerically 
resolved. In some cases, these simply saturate. In other cases, the 
nonlinear development is much more dramatic, leading to the 
formation of inward-moving density clumps. Cases with high 
relative field strengths (low plasma-/?) develop especially strong 
clumps. Such a clump grows exponentially with a growth time 
of a few orbits. The growth is tied to inhomogeneities in the 
strength of the outflow along different magnetic field lines. It 
is accompanied by a nose-shaped deformation of the magnetic 
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Figure 9. Snapshot of run J31l50 during the formation of the clump. Left panel: density (black/white) and field-parallel mass flux 
(green/pink colors). Right panel, azimuthal velocity (sub-Keplerian in blue and super-Keplerian in red) and azimuthal inclination of 
the magnetic field (brown/turquoise colors). The horizontal velocity, as measured in a frame moving with the mean velocity at the 
midplane, is visualized with small wedges whose areas are proportional to the velocity modulus. Selected magnetic field lines and 
the disk's scale height are overplotted in blue. See Fig. [5]: for a snapshot in the evolved state. 



^29,00 




Figure 8. Snapshots of the density in runs with fa = 1 and 
i - 45° for different values of tm- The second panel from the 
top shows the standard case (/3U45, compare Fig. [4]). The solid 
and dashed lines indicate the scale height and the sonic surface, 
respectively. 



field, with the clump located at the tip of the nose. Such a mas- 
sive instability has been proposed before in a simpler model of 




time [t =27i/n ] 



Figure 10. Growth of the clump in run /3U50. As a measure for 
the "dumpiness", we use the range of poloidal field inclinations 
at a fixed height in the disk. 



Lin et al. 



the disk-flow connection ( Agapitou 2000 ). Clumps may be a rea- 
son for the variability seen in observations of accreting systems, 
for instance the hard state noise in X-ray binaries (e.g 
|2000| ). 

At very low field strengths, the solutions become very differ- 
ent. "Classical" axisymmetric MRI modes develop if the field is 
vertical. The modes grow strong within roughly an orbit. With 
an inclined field, and despite reflective symmetry, the field tends 
to become extremely flat near the midplane. This then leads to 
magnetic reconnection. 

A major constraint of the simulations presented here is that 
3D effects are ignored. How do the 2D stripes in the wind 
turn into 3D structures? We circumvented the problem of how 
the material lost in the outflow is replenished. In reality, this 
is presumably achieved through 3D interchange processes. As 
suggested by 3D simulations of "magnet ically arrested" flows 
( Narayan et al.|2003|[Igumenshchev|2008| ), it is likely that high- 
amplitude clumps also form in the 3D case. This would also pro- 
vide justification for the "clumpy field" ac cretion of magnetic 
flux proposed by Spruit & Uzdensky (2005 ). 
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time [2n/fl D ] 

Figure 11. Evolution of various quantities along two different 
field lines (a) and (b) in run /3k50. From top to bottom: rotation 
velocity minus Keplerian value, mass flux, density, field incli- 
nation in azimuthal direction and field inclination in horizontal 
direction. 
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Figure 12. Snapshots of the azimuthal velocity (sub-Keplerian 
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lines (green) in two simulations with J3{ = 100: 1 = 90° in the 
top panel and 1 - 45° in the bottom panel. The initial azimuthal 
velocity was perturbed with white noise in both cases. 
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Figure 13. Snapshots of velocity, density and magnetic field 
lines in a simulation that begins in the evolved state of run /3k45 
(Fig. [4]) but with the magnetic field strength reduced by a factor 
10. Wedges depict the poloidal velocity in a frame moving with 
the mean horizontal velocity (v x ) = -3.8 /oAo- Red and blue 
colors represent super and sub-Keplerian azimuthal velocities, 
respectively. 
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